rm(list=ls())
knitr::opts_chunk$set(echo = TRUE, cache = TRUE)
library(tidyverse)
## ── Attaching packages ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.2.1 ✓ purrr 0.3.3
## ✓ tibble 2.1.3 ✓ dplyr 0.8.3
## ✓ tidyr 1.0.0 ✓ stringr 1.4.0
## ✓ readr 1.3.1 ✓ forcats 0.4.0
## ── Conflicts ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(deSolve)
library(rootSolve)
set.seed(13)
run_sims <- FALSE
The aim is to reproduce the model based results / figures.
Here is the original publication: https://www.nature.com/articles/s41467-017-00912-x
And here is the supplementary information, including the table of parameter values: https://static-content.springer.com/esm/art%3A10.1038%2Fs41467-017-00912-x/MediaObjects/41467_2017_912_MOESM1_ESM.pdf
Note that below I go with the notation from the table S1 and not the main text. However there is only one small difference between these: In the main text the notation is simplified to y_S_PB (which is y_SR_PB in Table S1) and y_S_SB (which is y_SO_SB in Table S1).
The published model is modified:
source("functions/bushplus_dynamic_model.r")
default_parameter_values <- c(
## maximum specific growth rates
g_max_CB = 0.05,
g_max_PB = 0.07,
g_max_SB = 0.1,
## half-saturation constants
k_PB_SR = 10,
k_SB_SO = 5,
k_CB_P = 0.2,
k_PB_P = 0.5,
k_SB_P = 0.5,
## half-inhibition constant
h_SR_CB = 300,
h_O_PB = 100,
h_O_SB = 100,
## yield (cells per micromole)
Y_SO_SB = 3.33e7,
Y_SR_PB = 1.25e7, ## Y_S_PB in the main text
y_P_CB = 1.67e8,
y_P_PB = 1.67e8,
y_P_SB = 1.67e8,
## oxygen production per microbial cell
P_CB = 6e-9,
## mortality rate
m_CB = 0.02,
m_PB = 0.028,
m_SB = 0.04,
## substrate diffusivity
a_S = 0.001,
a_O = 8e-4,
a_P = 0.01,
## background substrate concentration
back_SR = 300,
back_SO = 300,
back_O = 300,
back_P = 9.5,
## immigration rates (cells per time)
i_CB <- 0,
i_PB <- 0,
i_SB <- 0,
## oxidisation rate of reduced sulphur
c = 4e-5
)
Used to add noise and to ensure minimum abundance. Noise is added to the resource concentrations, and not to the organism abundances. Minimum organism abundance is set. I.e. organismal abundance cannot go below some very low level.
source("functions/default_event_definition.r")
This is when an event happens
default_event_interval <- 10
Default to no noise:
default_noise_sigma <- 0
Default to following minimums
default_minimum_abundances <- rep(1e2, 3)
names(default_minimum_abundances) <- c("CB", "PB", "SB")
default_sim_duration <- 2000
default_sim_sample_interval <- 1
Here we specify no forcing of a_0 (used unless set elsewhere).
default_log10a_series <- c(log10(default_parameter_values["a_O"]),
log10(default_parameter_values["a_O"]))
Initial conditions favouring anoxic
default_initial_state <- c(N_CB = 5e1,
N_PB = 1e7,
N_SB = 1e7,
SO = 300,
SR = 300,
O = 1e1,
P = 1e1)
source("functions/run_simulation.r")
source("functions/plot_dynamics.r")
source("functions/ss_expt_a_N.r")
source("functions/get_nonlinearity.r")
source("functions/get_hysteresis_range.r")
# dynamic_model = bushplus_dynamic_model
# event_definition = default_event_definition
# parameter_values = default_parameter_values
# event_interval = default_event_interval
# noise_sigma = default_noise_sigma
# minimum_abundances = default_minimum_abundances
# sim_duration = default_sim_duration
# sim_sample_interval = default_sim_sample_interval
# log10a_forcing = default_log10a_forcing
# initial_state = default_initial_state
# solver_method = "radau"
#
#
# times <- seq(0,
# sim_duration,
# by = sim_sample_interval)
#
# event_times <- seq(min(times),
# max(times),
# by = event_interval)
#
# log10a_forcing <- approxfun(x = log10a_forcing[,1],
# y = log10a_forcing[,2],
# method = "linear",
# rule = 2)
# assign("log10a_forcing", log10a_forcing, envir = .GlobalEnv)
#
# out <- as.data.frame(ode(y = initial_state,
# times = times,
# func = dynamic_model,
# parms = parameter_values,
# method = solver_method,
# events = list(func=event_definition,
# time=event_times)))
#
# sim_res <- run_simulation(dynamic_model = bushplus_dynamic_model,
# event_definition = default_event_definition,
# parameter_values = default_parameter_values,
# event_interval = default_event_interval,
# noise_sigma = default_noise_sigma,
# minimum_abundances = default_minimum_abundances,
# sim_duration = default_sim_duration,
# sim_sample_interval = default_sim_sample_interval,
# log10a_forcing = default_log10a_forcing,
# initial_state = default_initial_state,
# solver_method = "radau")
#sim_res <- run_simulation()
#plot_dynamics(sim_res)
Figure 2a, b in Bush et al.
init_state_favouring_anoxic <- c(N_CB = 5e1,
N_PB = 1e7,
N_SB = 1e7,
SO = 300,
SR = 300,
O = 1e1,
P = 1e1)
sim_res <- run_simulation(initial_state = init_state_favouring_anoxic)
plot_dynamics(simulation_result = sim_res)
#ggsave('bacteria_no_noise.png')
Population dynamics seem as in the paper. Resource dynamics fluctuate more than in the paper.
Figure 2a, b in Bush et al, but with noise
## add some noise
noise_sigma_0.01 <- 0.01
sim_res <- run_simulation(initial_state = init_state_favouring_anoxic,
noise_sigma = noise_sigma_0.01)
plot_dynamics(sim_res)
#ggsave("high_SB_PB.png")
Figure 2c, d in Bush et al.
init_state_favouring_oxic <- c(N_CB = 1e7,
N_PB = 1e5,
N_SB = 1e5,
SO = 500,
SR = 50,
O = 30,
P = 10)
sim_res <- run_simulation(initial_state = init_state_favouring_oxic)
plot_dynamics(sim_res)
##ggsave("high_CB.png")
Population and resource dynamics seen as those in the paper.
Figure 2c, d in Bush et al. but with noise.
sim_res <- run_simulation(initial_state = init_state_favouring_oxic,
noise_sigma = noise_sigma_0.01)
plot_dynamics(sim_res)
##ggsave("high_CB.png")
ssfind_minimum_abundances <- rep(0, 3)
names(ssfind_minimum_abundances) <- c("CB", "PB", "SB")
ssfind_simulation_duration <- 10000
ssfind_simulation_sampling_interval <- ssfind_simulation_duration
ssfind_event_interval <- ssfind_simulation_duration
ssfind_parameters <- default_parameter_values
# ssfind_init_state <- c(N_CB = NA,
# N_PB = NA,
# N_SB = NA,
# SO = 500,
# SR = 50,
# O = 30,
# P = 10)
##ssfind_init_state <- default_initial_state
ssfind_init_state <- init_state_favouring_anoxic
ssfind_init_state <- c(N_CB = NA, #this one is varied
N_PB = 1e8,
N_SB = 1e8,
SO = 200,
SR = 200,
O = 10,
P = 9.5)
grid_num_a <- 100
a_Os <- 10^seq(-4.5, -1.5, length=grid_num_a)
grid_num_N <- 2
initial_CBs <- 0
initial_PBs <- 0
initial_SBs <- 0
expt <- expand.grid(N_CB = initial_CBs,
N_PB = initial_PBs,
N_SB = initial_SBs,
a_O = a_Os)
ss_no_biology <- ss_by_a_N(expt)
saveRDS(ss_no_biology, "data/ss_no_biology.RDS")
ss_no_biology <- readRDS("data/ss_no_biology.RDS")
ss_no_biology %>%
select(a, N_CB, N_SB, N_PB) %>%
gather(key = Species, value=Quantity, 2:ncol(.)) %>%
ggplot() +
geom_path(aes(x = a, y = log10(Quantity), col = Species))
ss_no_biology %>%
select(a, SO, SR, O, P) %>%
gather(key = Species, value=Quantity, 2:ncol(.)) %>%
ggplot() +
geom_path(aes(x = a, y = log10(Quantity), col = Species))
initial_CBs <- 10^seq(0, 8, length=grid_num_N)
initial_PBs <- 0
initial_SBs <- 0
expt <- expand.grid(N_CB = initial_CBs,
N_PB = initial_PBs,
N_SB = initial_SBs,
a_O = a_Os)
ss_only_CB <- ss_by_a_N(expt)
saveRDS(ss_only_CB, "data/ss_only_CB.RDS")
ss_only_CB <- readRDS("data/ss_only_CB.RDS")
tempx1 <- ss_only_CB %>%
select(a, CB=N_CB, SB=N_SB, PB=N_PB) %>%
gather(key = Species, value=Quantity, 2:ncol(.))
tempx1 %>%
ggplot(aes(x = a, y = log10(Quantity), col = Species)) +
geom_path(size=0.8) +
ylab(expression(Log[10](Abundance))) +
xlab(expression(Log[10](Oxygen~diffusivity))) +
theme_light() +
ylim(0,10) +
geom_path(data = filter(tempx1, Species == "PB"),
aes(y=0.2), size = 0.8) +
theme(plot.title=element_text(size=10, hjust=0.5)) +
labs(title="Only CB")
ggsave("figures/CB.png", width = 3, height = 2)
ss_only_CB %>%
select(a, SO, SR, O, P) %>%
gather(key = Species, value=Quantity, 2:ncol(.)) %>%
ggplot() +
geom_path(aes(x = a, y = log10(Quantity), col = Species))
initial_CBs <- 0 # this one is varied
initial_PBs <- 0
initial_SBs <- 10^seq(0, 8, length=grid_num_N)
expt <- expand.grid(N_CB = initial_CBs,
N_PB = initial_PBs,
N_SB = initial_SBs,
a_O = a_Os)
ss_only_SB <- ss_by_a_N(expt)
saveRDS(ss_only_SB, "data/ss_only_SB.RDS")
ss_only_SB <- readRDS("data/ss_only_SB.RDS")
ss_only_SB %>%
select(a, N_CB, N_SB, N_PB) %>%
gather(key = Species, value=Quantity, 2:ncol(.)) %>%
ggplot() +
geom_path(aes(x = a, y = log10(Quantity), col = Species))
ss_only_SB %>%
select(a, SO, SR, O, P) %>%
gather(key = Species, value=Quantity, 2:ncol(.)) %>%
ggplot() +
geom_path(aes(x = a, y = log10(Quantity), col = Species))
initial_CBs <- 0 # this one is varied
initial_PBs <- 10^seq(0, 8, length=grid_num_N)
initial_SBs <- 0
expt <- expand.grid(N_CB = initial_CBs,
N_PB = initial_PBs,
N_SB = initial_SBs,
a_O = a_Os)
ss_only_PB <- ss_by_a_N(expt)
saveRDS(ss_only_PB, "data/ss_only_PB.RDS")
ss_only_PB <- readRDS("data/ss_only_PB.RDS")
tempx1 <- ss_only_PB %>%
select(a, CB=N_CB, SB=N_SB, PB=N_PB) %>%
gather(key = Species, value=Quantity, 2:ncol(.))
tempx1 %>%
ggplot(aes(x = a, y = log10(Quantity), col = Species)) +
geom_path(size = 0.8) +
theme_light() +
ylim(0,10) +
ylab(expression(Log[10](Abundance))) +
xlab(expression(Log[10](Oxygen~diffusivity))) +
geom_path(data = filter(tempx1, Species == "CB"),
aes(y=0.2), size = 0.8)+
geom_path(data = filter(tempx1, Species == "PB"),
size = 0.8) +
theme(plot.title=element_text(size=10, hjust=0.5)) +
labs(title="Only PB")
ggsave("figures/PB.png", width = 3, height = 2)
ss_only_PB %>%
select(a, SO, SR, O, P) %>%
gather(key = Species, value=Quantity, 2:ncol(.)) %>%
ggplot() +
geom_path(aes(x = a, y = log10(Quantity), col = Species))
initial_CBs <- 0 # this one is varied
initial_PBs <- 10^seq(0, 8, length=grid_num_N)
initial_SBs <- 0
expt <- expand.grid(N_CB = initial_CBs,
N_PB = initial_PBs,
N_SB = initial_SBs,
a_O = a_Os)
expt$N_SB <- expt$N_PB
ss_only_SB_PB <- ss_by_a_N(expt)
saveRDS(ss_only_SB_PB, "data/ss_only_SB_PB.RDS")
ss_only_SB_PB <- readRDS("data/ss_only_SB_PB.RDS")
ss_only_SB_PB %>%
select(a, N_CB, N_SB, N_PB) %>%
gather(key = Species, value=Quantity, 2:ncol(.)) %>%
ggplot() +
geom_path(aes(x = a, y = log10(Quantity), col = Species))
ss_only_SB_PB %>%
select(a, SO, SR, O, P) %>%
gather(key = Species, value=Quantity, 2:ncol(.)) %>%
ggplot() +
geom_path(aes(x = a, y = log10(Quantity), col = Species))
initial_CBs <- 10^seq(0, 8, length=grid_num_N)
initial_PBs <- 0
initial_SBs <- 0
expt <- expand.grid(N_CB = initial_CBs,
N_PB = initial_PBs,
N_SB = initial_SBs,
a_O = a_Os)
expt$N_SB <- expt$N_CB
ss_only_SB_CB <- ss_by_a_N(expt)
saveRDS(ss_only_SB_CB, "data/ss_only_SB_CB.RDS")
ss_only_SB_CB <- readRDS("data/ss_only_SB_CB.RDS")
ss_only_SB_CB %>%
select(a, N_CB, N_SB, N_PB) %>%
gather(key = Species, value=Quantity, 2:ncol(.)) %>%
ggplot() +
geom_path(aes(x = a, y = log10(Quantity), col = Species))
ss_only_SB_CB %>%
select(a, SO, SR, O, P) %>%
gather(key = Species, value=Quantity, 2:5) %>%
ggplot() +
geom_path(aes(x = a, y = log10(Quantity), col = Species))
initial_CBs <- 0 # this one is varied
initial_PBs <- 10^seq(0, 8, length=grid_num_N)
initial_SBs <- 0
expt <- expand.grid(N_CB = initial_CBs,
N_PB = initial_PBs,
N_SB = initial_SBs,
a_O = a_Os)
expt$N_CB <- expt$N_PB
ss_only_CB_PB <- ss_by_a_N(expt)
saveRDS(ss_only_CB_PB, "data/ss_only_CB_PB.RDS")
ss_only_CB_PB <- readRDS("data/ss_only_CB_PB.RDS")
ss_only_CB_PB %>%
select(a, N_CB, N_SB, N_PB) %>%
gather(key = Species, value=Quantity, 2:ncol(.)) %>%
ggplot() +
geom_path(aes(x = a, y = log10(Quantity), col = Species))
ss_only_SB_PB %>%
select(a, SO, SR, O, P) %>%
gather(key = Species, value=Quantity, 2:ncol(.)) %>%
ggplot() +
geom_path(aes(x = a, y = log10(Quantity), col = Species))
initial_CBs <- 10^seq(0, 8, length=grid_num_N)
initial_PBs <- 1e5
initial_SBs <- 1e5
expt <- expand.grid(N_CB = initial_CBs,
N_PB = initial_PBs,
N_SB = initial_SBs,
a_O = a_Os)
ss_CB_SB_PB <- ss_by_a_N(expt)
saveRDS(ss_CB_SB_PB, "data/ss_CB_SB_PB.RDS")
ss_CB_SB_PB <- readRDS("data/ss_CB_SB_PB.RDS")
ss_CB_SB_PB %>%
select(a, CB=N_CB, SB=N_SB, PB=N_PB) %>%
gather(key = Species, value=Quantity, 2:ncol(.)) %>%
ggplot() +
geom_path(aes(x = a, y = log10(Quantity), col = Species),
size=0.8) +
ylab(expression(Log[10](Abundance))) +
xlab(expression(Log[10](Oxygen~diffusivity))) +
theme_light() +
ylim(0,10) +
theme(plot.title=element_text(size=10, hjust=0.5)) +
labs(title="CB, SB, & PB")
ggsave("figures/CB_SB_PB.png", width = 3, height = 2)
ss_CB_SB_PB %>%
select(a, SO, SR, O, P) %>%
gather(key = Species, value=Quantity, 2:ncol(.)) %>%
ggplot() +
geom_path(aes(x = a, y = log10(Quantity), col = Species))
all_ss <- bind_cols(readRDS("data/ss_no_biology.RDS"),
composition = rep("No organisms", grid_num_a*grid_num_N),
direction = c(rep("up", grid_num_a),
rep("down", grid_num_a)))
all_ss <- bind_rows(all_ss,
bind_cols(readRDS("data/ss_only_CB.RDS"),
composition = rep("CB", grid_num_a*grid_num_N),
direction = c(rep("up", grid_num_a),
rep("down", grid_num_a))))
all_ss <- bind_rows(all_ss,
bind_cols(readRDS("data/ss_only_SB.RDS"),
composition = rep("SB", grid_num_a*grid_num_N),
direction = c(rep("up", grid_num_a),
rep("down", grid_num_a))))
all_ss <- bind_rows(all_ss,
bind_cols(readRDS("data/ss_only_PB.RDS"),
composition = rep("PB", grid_num_a*grid_num_N),
direction = c(rep("up", grid_num_a),
rep("down", grid_num_a))))
all_ss <- bind_rows(all_ss,
bind_cols(readRDS("data/ss_only_SB_PB.RDS"),
composition = rep("SB-PB", grid_num_a*grid_num_N),
direction = c(rep("up", grid_num_a),
rep("down", grid_num_a))))
all_ss <- bind_rows(all_ss,
bind_cols(readRDS("data/ss_only_SB_CB.RDS"),
composition = rep("SB-CB", grid_num_a*grid_num_N),
direction = c(rep("up", grid_num_a),
rep("down", grid_num_a))))
all_ss <- bind_rows(all_ss,
bind_cols(readRDS("data/ss_only_CB_PB.RDS"),
composition = rep("CB-PB", grid_num_a*grid_num_N),
direction = c(rep("up", grid_num_a),
rep("down", grid_num_a))))
all_ss <- bind_rows(all_ss,
bind_cols(readRDS("data/ss_CB_SB_PB.RDS"),
composition = rep("CB-SB-PB",
grid_num_a*grid_num_N),
direction = c(rep("up", grid_num_a),
rep("down", grid_num_a))))
x_long <- all_ss %>%
gather(key = Species, value = Quantity, 2:8)
x_wide <- x_long %>%
select(-initial_N_CB) %>%
spread(key = direction, value=Quantity, drop=T)
resp_summary <- x_wide %>%
group_by(composition, Species) %>%
summarise(hyst_1 = mean(abs(log10(up) - log10(down))),
hyst_2 = get_hysteresis_range(up, down, a),
nl_up = get_nonlinearity(a, log10(up)),
nl_down = get_nonlinearity(a, log10(down))) %>%
ungroup() %>%
gather(key = resp_measure, value = value, 3:6) %>%
mutate(composition = fct_relevel(composition,
"No organisms",
"CB",
"SB",
"PB",
"CB-PB",
"SB-CB",
"SB-PB",
"CB-SB-PB"),
num_func_groups = case_when(composition =="No organisms" ~ 0,
composition == "CB" ~ 1,
composition == "SB" ~ 1,
composition == "PB" ~ 1,
composition == "CB-PB" ~ 2,
composition == "SB-CB" ~ 2,
composition == "SB-PB" ~ 2,
composition == "CB-SB-PB" ~ 3
),
var_type = ifelse(str_sub(Species, 1, 1)=="N", "Organism", "Substrate")
)
saveRDS(resp_summary, "data/resp_summary.RDS")
resp_summary <- readRDS("data/resp_summary.RDS")
resp_summary %>%
ggplot() +
geom_point(aes(x = composition, y = value)) +
facet_grid(Species ~ resp_measure, scales = "free") +
coord_flip()
resp_summ1 <- resp_summary %>%
filter(resp_measure %in% c("nl_down", "nl_up")) %>%
group_by(composition, Species, num_func_groups, var_type) %>%
summarise(ave_nl = mean(value))
resp_summ1 %>%
ggplot() +
geom_point(aes(composition, ave_nl)) +
facet_wrap( ~ Species)
focus_on <- "Substrate"
ave_nl_substrate <- resp_summ1 %>%
filter(var_type == focus_on) %>%
group_by(composition) %>%
summarise(ave_ave_nl = mean(ave_nl))
ave_nl_substrate %>%
ggplot() +
geom_point(aes(composition, ave_ave_nl)) +
ggtitle(paste("Nonlinearity in", focus_on, "response."))
focus_on <- "Organism"
ave_nl_organism <- resp_summ1 %>%
filter(var_type == focus_on) %>%
group_by(composition) %>%
summarise(ave_ave_nl = sum(ave_nl)/unique(num_func_groups))
ave_nl_organism %>%
ggplot() +
geom_point(aes(composition, ave_ave_nl)) +
ggtitle(paste("Nonlinearity in", focus_on, "response."))
## Warning: Removed 1 rows containing missing values (geom_point).
focus_on <- "System"
ave_nl_system <- resp_summ1 %>%
group_by(composition) %>%
summarise(ave_ave_nl = sum(ave_nl) / (unique(num_func_groups)+4))
ave_nl_system %>%
ggplot() +
geom_point(aes(composition, ave_ave_nl)) +
ggtitle(paste("Nonlinearity in", focus_on, "response."))
ggsave("figures/system_nonlinearity.png", width = 3, height = 2)
focus_on <- "Substrate"
resp_summ3 <- resp_summary %>%
filter(resp_measure %in% c("hyst_1", "hyst_2"),
var_type == focus_on) %>%
group_by(composition, num_func_groups, resp_measure, var_type) %>%
summarise(ave_hyst = mean(value))
resp_summ3 %>%
ggplot() +
geom_point(aes(x = composition, y = ave_hyst)) +
facet_wrap( ~ resp_measure, scales="free") +
ggtitle(paste("Hysteresis in", focus_on, "response."))
focus_on <- "Organism"
resp_summ3 <- resp_summary %>%
filter(resp_measure %in% c("hyst_1", "hyst_2"),
var_type == focus_on) %>%
group_by(composition, num_func_groups, resp_measure, var_type) %>%
summarise(ave_hyst = sum(value) / unique(num_func_groups))
resp_summ3 %>%
ggplot() +
geom_point(aes(x = composition, y = ave_hyst)) +
facet_wrap( ~ resp_measure, scales="free") +
ggtitle(paste("Hysteresis in", focus_on, "response."))
## Warning: Removed 2 rows containing missing values (geom_point).
resp_summ4 <- resp_summary %>%
filter(resp_measure %in% c("hyst_1", "hyst_2")) %>%
group_by(resp_measure, composition) %>%
summarise(mean_hyst = sum(value) / (unique(num_func_groups)+4))
resp_summ4 %>%
ggplot() +
geom_point(aes(x = composition, y = mean_hyst)) +
facet_wrap( ~ resp_measure) +
ggtitle("Hysteresis in system response.")
resp_summ5 <- resp_summ4 %>%
group_by(composition) %>%
summarise(mean_mean_hyst = mean(mean_hyst))
resp_summ5 %>%
ggplot() +
geom_point(aes(x = composition, y = mean_mean_hyst)) +
xlab("Functional groups present") +
ylab("Amount of hysteresis") +
ylim(0,2) +
theme_light() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
ggsave("figures/composition_system_hysteresis.png", width = 2.5, height = 2.5)
Three parameters set the strength of inhibition (half-inhibition constants):
These work by the inhibition function:
\(\frac{1}{1 + \frac{X}{h_X}}\)
inhibition
## function (X, h_X)
## 1/(1 + (X/h_X))
Plot some examples:
default_parameter_values["h_SR_CB"]
## h_SR_CB
## 300
SR <- 10^seq(1.4, 2.5, length = 100)
ggplot() +
geom_line(aes(x = SR, y = inhibition(SR, default_parameter_values["h_SR_CB"]))) +
ylab("Multiplicative effect\non CB growth rate")
default_parameter_values["h_O_PB"]
## h_O_PB
## 100
O <- 10^seq(-0.2, 2.5, length = 100)
ggplot() +
geom_line(aes(x = O, y = inhibition(O, default_parameter_values["h_O_PB"]))) +
ylab("Multiplicative effect\non PB growth rate")
default_parameter_values["h_O_SB"]
## h_O_SB
## 100
O <- 10^seq(-0.2, 2.5, length = 100)
ggplot() +
geom_line(aes(x = O, y = inhibition(O, default_parameter_values["h_O_SB"]))) +
ylab("Multiplicative effect\non SB growth rate")
proportion_change <- 2
default_parameter_values["h_O_SB"]
## h_O_SB
## 100
O <- 10^seq(-0.2, 2.5, length = 100)
ggplot() +
geom_line(aes(x = O, y = inhibition(O, proportion_change*default_parameter_values["h_O_SB"]))) +
ylab("Multiplicative effect\non SB growth rate")
Increasing the inhibition half saturation constants weakens the inhibition.
ssfind_minimum_abundances <- default_minimum_abundances
ssfind_simulation_duration <- 10000
ssfind_simulation_sampling_interval <- ssfind_simulation_duration
ssfind_event_interval <- ssfind_simulation_duration
ssfind_parameters <- default_parameter_values
ssfind_init_state <- c(N_CB = NA, #this one is varied
N_PB = NA,
N_SB = NA,
SO = 200,
SR = 200,
O = 10,
P = 9.5)
grid_num_mult <- 20
grid_num_N <- 2
grid_num_a <- 100
a_Os <- 10^seq(-4.5, -1.5, length=grid_num_a)
inhib_const_mults <- 2^seq(-1, 0.4, length=grid_num_mult)
initial_CBs <- 10^seq(0, 8, length=grid_num_N)
initial_PBs <- 1e5
initial_SBs <- 1e5
expt <- expand.grid(N_CB = initial_CBs,
N_PB = initial_PBs,
N_SB = initial_SBs,
inhib_const_mult = inhib_const_mults,
a_O = a_Os)
ss_h1 <- ss_by_a_h_N(expt)
saveRDS(ss_h1, "data/ss_h1.RDS")
ss_h1 <- readRDS("data/ss_h1.RDS")
inhib_const_mults
## [1] 0.5000000 0.5262004 0.5537737 0.5827919 0.6133306 0.6454696 0.6792927
## [8] 0.7148882 0.7523489 0.7917725 0.8332620 0.8769256 0.9228772 0.9712366
## [15] 1.0221302 1.0756906 1.1320576 1.1913783 1.2538074 1.3195079
focal_inhibition_multiplier <- inhib_const_mults[14]
ss_h1 %>%
filter(inhib_const_mult == focal_inhibition_multiplier) %>%
select(a, N_CB, N_SB, N_PB) %>%
gather(key = Species, value=Quantity, 2:ncol(.)) %>%
ggplot() +
geom_path(aes(x = a, y = log10(Quantity), col = Species))
ss_h1 %>%
filter(inhib_const_mult == focal_inhibition_multiplier) %>%
select(a, SO, SR, O, P) %>%
gather(key = Species, value=Quantity, 2:ncol(.)) %>%
ggplot() +
geom_path(aes(x = a, y = log10(Quantity), col = Species))
## YES! 0.5 does mean twice as strong inhibition
focal_inhibition_multiplier <- 0.5
ss_h1 %>%
filter(inhib_const_mult == focal_inhibition_multiplier) %>%
select(a, N_CB, N_SB, N_PB) %>%
gather(key = Species, value=Quantity, 2:ncol(.)) %>%
ggplot() +
geom_path(aes(x = a, y = log10(Quantity), col = Species))
ss_h1 %>%
filter(inhib_const_mult == focal_inhibition_multiplier) %>%
select(a, SO, SR, O, P) %>%
gather(key = Species, value=Quantity, 2:ncol(.)) %>%
ggplot() +
geom_path(aes(x = a, y = log10(Quantity), col = Species))
## YES! 2 does mean half as strong inhibition
focal_inhibition_multiplier <- 2^0.4
ss_h1 %>%
filter(inhib_const_mult == focal_inhibition_multiplier) %>%
select(a, N_CB, N_SB, N_PB) %>%
gather(key = Species, value=Quantity, 2:ncol(.)) %>%
ggplot() +
geom_path(aes(x = a, y = log10(Quantity), col = Species))
ss_h1 %>%
filter(inhib_const_mult == focal_inhibition_multiplier) %>%
select(a, SO, SR, O, P) %>%
gather(key = Species, value=Quantity, 2:ncol(.)) %>%
ggplot() +
geom_path(aes(x = a, y = log10(Quantity), col = Species))
x_long <- ss_h1 %>%
mutate(direction = ifelse(initial_N_CB==1, "up", "down")) %>%
gather(key = Species, value = Quantity, 3:9)
x_wide <- x_long %>%
select(-initial_N_CB) %>%
spread(key = direction, value=Quantity, drop=T)
resp_summary <- x_wide %>%
group_by(inhib_const_mult, Species) %>%
summarise(hyst_1 = mean(abs(log10(up) - log10(down))),
hyst_2 = get_hysteresis_range(up, down, a),
nl_up = get_nonlinearity(a, log10(up)),
nl_down = get_nonlinearity(a, log10(down))) %>%
ungroup() %>%
gather(key = resp_measure, value = value, 3:6) %>%
mutate(
var_type = ifelse(str_sub(Species, 1, 1)=="N", "Organism", "Substrate")
)
resp_summ1 <- resp_summary %>%
filter(resp_measure %in% c("nl_down", "nl_up")) %>%
group_by(inhib_const_mult, Species, var_type) %>%
summarise(ave_nl = mean(value))
resp_summ1 %>%
ggplot() +
geom_point(aes(inhib_const_mult, ave_nl)) +
facet_wrap( ~ Species)
focus_on <- "Substrate"
ave_nl_substrate <- resp_summ1 %>%
filter(var_type == focus_on) %>%
group_by(inhib_const_mult) %>%
summarise(ave_ave_nl = mean(ave_nl))
ave_nl_substrate %>%
ggplot() +
geom_point(aes(inhib_const_mult, ave_ave_nl)) +
ggtitle(paste("Nonlinearity in", focus_on, "response."))
focus_on <- "Organism"
ave_nl_organism <- resp_summ1 %>%
filter(var_type == focus_on) %>%
group_by(inhib_const_mult) %>%
summarise(ave_ave_nl = mean(ave_nl))
ave_nl_organism %>%
ggplot() +
geom_point(aes(inhib_const_mult, ave_ave_nl)) +
ggtitle(paste("Nonlinearity in", focus_on, "response."))
focus_on <- "System"
ave_nl_system <- resp_summ1 %>%
group_by(inhib_const_mult) %>%
summarise(ave_ave_nl = mean(ave_nl))
ave_nl_system %>%
ggplot() +
geom_point(aes(inhib_const_mult, ave_ave_nl)) +
ggtitle(paste("Nonlinearity in", focus_on, "response."))
ggsave("figures/system_nonlinearity.png", width = 3, height = 2)
focus_on <- "Substrate"
resp_summ3 <- resp_summary %>%
filter(resp_measure %in% c("hyst_1", "hyst_2"),
var_type == focus_on) %>%
group_by(inhib_const_mult, resp_measure, var_type) %>%
summarise(ave_hyst = mean(value))
resp_summ3 %>%
ggplot() +
geom_point(aes(x = inhib_const_mult, y = ave_hyst)) +
facet_wrap( ~ resp_measure, scales="free") +
ggtitle(paste("Hysteresis in", focus_on, "response."))
focus_on <- "Organism"
resp_summ3 <- resp_summary %>%
filter(resp_measure %in% c("hyst_1", "hyst_2"),
var_type == focus_on) %>%
group_by(inhib_const_mult, resp_measure, var_type) %>%
summarise(ave_hyst = mean(value))
resp_summ3 %>%
ggplot() +
geom_point(aes(x = inhib_const_mult, y = ave_hyst)) +
facet_wrap( ~ resp_measure, scales="free") +
ggtitle(paste("Hysteresis in", focus_on, "response."))
resp_summ4 <- resp_summary %>%
filter(resp_measure %in% c("hyst_1", "hyst_2")) %>%
group_by(resp_measure, inhib_const_mult) %>%
summarise(mean_hyst = mean(value))
resp_summ4 %>%
ggplot() +
geom_point(aes(x = inhib_const_mult, y = mean_hyst)) +
facet_wrap( ~ resp_measure) +
ggtitle("Hysteresis in system response.")
resp_summ5 <- resp_summ4 %>%
group_by(inhib_const_mult) %>%
summarise(mean_mean_hyst = mean(mean_hyst))
resp_summ5 %>%
ggplot() +
geom_point(aes(x = inhib_const_mult, y = mean_mean_hyst)) +
xlab("Organismal tolerance") +
ylab("Amount of hysteresis") +
theme_light() +
ylim(0,2)
ggsave("figures/inhibition_system_hysteresis.png", width = 2.5, height = 2)
## anoxic initial state
init_state <- init_state_favouring_anoxic
## simulation timing
simulation_duration <- 5000
simulation_sampling_interval <- 10
event_interval_t <- 100
## low-high-low oxgygen diffusivity
lhl_log10a_series <- c(-8, -8, -2, -2, -2, -2, -2, -8, -8)
## run simulation
sim_res <- run_simulation(initial_state = init_state_favouring_anoxic,
log10a_series = lhl_log10a_series,
sim_duration = simulation_duration,
sim_sample_interval = simulation_sampling_interval,
event_interval = event_interval_t)
saveRDS(sim_res, "data/lhl_with_defaults.RDS")
lhl_with_defaults <- readRDS("data/lhl_with_defaults.RDS")
plot_dynamics(lhl_with_defaults)
lhl_with_defaults$result %>%
filter(between(time, 1000, 73000)) %>%
ggplot() +
geom_path(aes(x=a, y=log10(O)),
arrow = arrow(type = "open", angle = 30, length = unit(0.1, "inches")))
## anoxic initial state
init_state <- init_state_favouring_oxic
## simulation timing
simulation_duration <- 50000
simulation_sampling_interval <- 10
event_interval_t <- 100
## low-high-low oxgygen diffusivity
hlh_log10a_series <- c(-2, -2, -10, -14, -18, -22, -2, -2, -2)
## run simulation
sim_res <- run_simulation(initial_state = init_state_favouring_oxic,
log10a_series = hlh_log10a_series,
sim_duration = simulation_duration,
sim_sample_interval = simulation_sampling_interval,
event_interval = event_interval_t)
saveRDS(sim_res, "data/hlh_with_defaults.RDS")
hlh_with_defaults <- readRDS("data/hlh_with_defaults.RDS")
plot_dynamics(hlh_with_defaults)
hlh_with_defaults$result %>%
filter(between(time, 1000, 100000)) %>%
ggplot() +
geom_path(aes(x=a, y=log10(O)),
arrow = arrow(type = "open", angle = 30, length = unit(0.1, "inches")))
init_state_anoxic_2 <- c(N_CB = 1e1,
N_PB = 1e8,
N_SB = 1e8,
SO = 300,
SR = 300,
O = 1e-5,
P = 1e1)
RS <- runsteady(y = init_state_anoxic_2,
fun = bushplus_dynamic_model,
parms = default_parameter_values,
times = c(0, 1e5),
stol=1)
RS$y
RS <- runsteady(y = init_state_favouring_oxic,
fun = bushplus_dynamic_model,
parms = default_parameter_values,
times = c(0, 1e5),
stol=1)
RS$y
Note that this section not working / updated since the code was made into functions
noise_sigma <- 0.0
get_final_states_N_CB <- function(x) {
#give N_CB, N_SB, N_PB and a_O values and put these into state and parameters
state["N_CB"] <- x["N_CB"]
parameters["a_O"] <- x["a_O"]
## update forcing function
a_forcing1 <- matrix(ncol=2, byrow=T, data=c(0,
log10(parameters["a_O"]),
max(times),
log10(parameters["a_O"])))
log10a_forcing2 <- approxfun(x = a_forcing1[,1],
y = a_forcing1[,2],
method = "linear",
rule = 2)
## assign the changed a_forcing2 into the global environment, from where the model gets a_forcing2
assign("log10a_forcing2", log10a_forcing2, envir = .GlobalEnv)
as.data.frame(ode(y = state,
times = times,
func = micro_mod_var,
parms = parameters,
method = "radau",
events = list(func=noise_events, time=times)))[2,-1]
#gives you values for the last time step
}
times <- seq(0, 50000, by = 50000)
state <- c(N_CB = NA, #this one is varied
N_PB = 1e8,
N_SB = 1e8,
SO = 200,
SR = 200,
O = 10,
P = 9.5)
#different combinations of N_CB and oxygen diffusability
grid_num <- 50
expt <- expand.grid(N_CB = 10^seq(0, 8, length=grid_num),
a_O = 10^seq(-3.5, -2, length=grid_num))
result <- apply(expt, 1, function(x) get_final_states_N_CB(x)) %>%
tibble() %>%
unnest() %>%
mutate(initial_N_CB=expt$N_CB,
a_O=expt$a_O,
highlow=ifelse(N_CB>2e8, 'high', 'low'))
result <- result %>%
mutate(N_CB = ifelse(N_CB < minimum_abundance, minimum_abundance, N_CB))
result$highlow <- as.factor(result$highlow)
saveRDS(result, "data/CB_oxdiff_hyst.RDS")
result <- readRDS("data/CB_oxdiff_hyst.RDS")
#plot N_CB / P_CB against a_0
p <- ggplot(result) +
geom_line(aes(x=log10(a_O), y=log10(N_CB), color=highlow), size=1.25) +
ggtitle("Stable states of the system")+
xlab('log[oxygen diffusivity]') +
ylab('Initial cyanobacteria abundance')
p <- p + geom_point(data=result, aes(x=log10(a_O), y=log10(initial_N_CB), color=highlow)) + labs(color = "Final CB concentration")
plot(p)
#ggsave("initial conditions_1.png")
#show hysterisis effects for all the parameters
#plot oxygen against diffusability
ggplot(result) +
geom_point(aes(x=log10(a_O), y=(log10(O))))+
xlab('log[oxygen diffusivity]') +
ylab('log[oxygen concentration]')+
geom_segment(aes(x=log10(result$a_O[839]), y=log10(result$O[839]), xend=log10(result$a_O[838]), yend=log10(result$O[838])), arrow=arrow(), colour="blue", size=0.5)+
geom_segment(aes(x=log10(result$a_O[1761]), y=log10(result$O[1761]), xend=log10(result$a_O[1762]), yend=log10(result$O[1762])), arrow=arrow(), colour="blue", size=0.5)+
geom_segment(aes(x=-4.25, y=-0.45, xend=-3.75, yend=0.1), arrow=arrow(), colour="blue", size=0.5)+
geom_curve(aes(x = -3.75, y = 1.9, xend = -4.25, yend = 1.5), arrow=arrow(), colour='blue', size = 0.5, curvature = -0.2)
#ggsave("hysterisis_CB_1.png")
ggplot(result) +
geom_point(aes(x=log10(a_O), y=(log10(P))))+
ggtitle("Hysterisis effect of phoshate concentration")+
xlab('log[oxygen diffusivity]') +
ylab('log[phosphate concentration]')
#ggsave("hysterisis_CB_2.png")
ggplot(result) +
geom_point(aes(x=log10(a_O), y=(log10(SO))))+
ggtitle("Hysterisis effect of oxidized sulfur concentration")+
xlab('log[oxygen diffusivity]') +
ylab('log[oxydized sulfur concentration]')
#ggsave("hysterisis_CB_3.png")
ggplot(result) +
geom_point(aes(x=log10(a_O), y=(log10(SR))))+
ggtitle("Hysterisis effect of reduced sulfur concentration")+
xlab('log[oxygen diffusivity]') +
ylab('log[reduced sulfur] concentration]')
#ggsave("hysterisis_CB_4.png")
These graphs show an at least qualitative match to the result in the paper.
noise_sigma <- 0.00
get_final_states_N_PB <- function(x) {
#give in N_CB, N_SB, N_PB and a_O values and put these into state and parameters
state["N_PB"] <- x["var"]
state["N_SB"] <- x["var"]
parameters["a_O"] <- x["a_O"]
## update forcing function
a_forcing1 <- matrix(ncol=2, byrow=T, data=c(0,
log10(parameters["a_O"]),
max(times),
log10(parameters["a_O"])))
log10a_forcing2 <- approxfun(x = a_forcing1[,1],
y = a_forcing1[,2],
method = "linear",
rule = 2)
## assign the changed a_forcing2 into the global environment, from where the model gets a_forcing2
assign("log10a_forcing2", log10a_forcing2, envir = .GlobalEnv)
as.data.frame(ode(y = state,
times = times,
func = micro_mod_var,
parms = parameters,
method = "radau", events = list(func=noise_events, time=times)))[2,-1]
#gives you values for the last time step
}
times <- seq(0, 50000, by = 5000)
state <- c(N_CB = 5,
N_PB = NA,#this one is varied
N_SB = NA,
SO = 200,
SR = 200,
O = 10,
P = 9.5)
#different combinations of N_PB and oxygen diffusivity
grid_num <- 50
expt <- expand.grid(var=10^seq(0, 7, length=grid_num),
a_O = 10^seq(-3, -2, length=grid_num))
result2 <- apply(expt, 1, function(x) get_final_states_N_PB(x)) %>%
tibble() %>%
unnest() %>%
mutate(initial_N_PB=expt$var,
initial_N_SB=expt$var,
a_O=expt$a_O,
highlow=ifelse(N_PB>1e8, 'high', 'low'))
result2 <- result2 %>%
mutate(N_PB = ifelse(N_PB < minimum_abundance, minimum_abundance, N_PB))
result2$highlow <- as.factor(result2$highlow)
saveRDS(result2, "data/PBSB_oxdiff_hyst.RDS")
result2 <- readRDS("data/PBSB_oxdiff_hyst.RDS")
#plot N_PB against a_0
p <- ggplot(result2) +
geom_line(aes(x=log10(a_O), y=log10(N_PB), color=highlow), size=1.25) +
ggtitle("Stable states of the system")+
xlab('log[oxygen diffusivity]') +
ylab('Initial PS and SR bacteria concentration')
p <- p + geom_point(data=result2,
aes(x=log10(a_O), y=log10(initial_N_PB), color=highlow))
p <- p + labs(color = "Final PB and SB concentration")
plot(p)
#ggsave("initial conditions_2.png")
Evaluation of the switch: Get the functions:
source('C:/Users/probs/Desktop/Sem 6/Internship/functions_analysis_EWS_2.R')
Get state shifts and calculate parameters:
x3 <- 30 #size of wished rolling window (percentage btw 0 and 100)
x5 <- 50 #size of wished rolling window (percentage btw 0 and 100)
#first switch
name_1='anoxic_oxic 1'
trans1 <- data.frame(out_t %>% filter(time>87500 & time<92500)) #for analysis
trans1_f <- data.frame(out_f %>% filter(time>87500 & time<92500)) #for plotting
n_13 <- round(length(trans1$time)*x3/100) #size of rolling window
t_13 <- length(trans1$time) - n_13 + 1 #number of rolling windows
n_15 <- round(length(trans1$time)*x5/100) #size of rolling window
t_15 <- length(trans1$time) - n_15 + 1 #number of rolling windows
out_df_1_N_CB <- data.frame(trans1$N_CB)
colnames(out_df_1_N_CB) <- c("N_CB")
out_df_1_N_PB <- data.frame(trans1$N_PB)
colnames(out_df_1_N_PB) <- c("N_PB")
out_df_1_N_SB <- data.frame(trans1$N_SB)
colnames(out_df_1_N_SB) <- c("N_SB")
ggplot(aes(x=time, y=trans_quantity, col=species), data=trans1_f) +
ggtitle("Anoxic_oxic 1")+
xlab('Time') +
ylab('Trans_quantity') +
facet_wrap(~var_type, scales = "free", ncol = 1) + #create one plot for organisms and one for substrate
geom_line()+
labs(color = "Organism/Substrate")
ggsave("anoxic_oxic 1.png")
#second switch
name_2='oxic_anoxic 1'
trans2 <- data.frame(out_t %>% filter(time>335000 & time<360000))
trans2_f <- data.frame(out_f %>% filter(time>335000 & time<360000))
n_23 <- round(length(trans2$time)*x3/100)
t_23 <- length(trans2$time) - n_23 + 1
n_25 <- round(length(trans2$time)*x5/100)
t_25 <- length(trans2$time) - n_25 + 1
out_df_2_N_CB <- data.frame(trans2$N_CB)
colnames(out_df_2_N_CB) <- c("N_CB")
out_df_2_N_PB <- data.frame(trans2$N_PB)
colnames(out_df_2_N_PB) <- c("N_PB")
out_df_2_N_SB <- data.frame(trans2$N_SB)
colnames(out_df_2_N_SB) <- c("N_SB")
ggplot(aes(x=time, y=trans_quantity, col=species), data=trans2_f) +
ggtitle("Oxic_anoxic 1")+
xlab('Time') +
ylab('Trans_quantity') +
facet_wrap(~var_type, scales = "free", ncol = 1) + #create one plot for organisms and one for substrate
geom_line()+
labs(color = "Organism/Substrate")
ggsave("oxic_anoxic 1.png")
#third switch
name_3='oxic_anoxic 2'
trans3 <- data.frame(out2_t %>% filter(time>165000 & time<185000))
trans3_f <- data.frame(out2_f %>% filter(time>165000 & time<185000))
n_33 <- round(length(trans3$time)*x3/100)
t_33 <- length(trans3$time) - n_33 + 1
n_35 <- round(length(trans3$time)*x5/100)
t_35 <- length(trans3$time) - n_35 + 1
out2_df_1_N_CB <- data.frame(trans3$N_CB)
colnames(out2_df_1_N_CB) <- c("N_CB")
out2_df_1_N_SB <- data.frame(trans3$N_SB)
colnames(out2_df_1_N_SB) <- c("N_SB")
out2_df_1_N_PB <- data.frame(trans3$N_PB)
colnames(out2_df_1_N_PB) <- c("N_PB")
ggplot(aes(x=time, y=trans_quantity, col=species), data=trans3_f) +
ggtitle("Oxic_anoxic 2")+
xlab('Time') +
ylab('Trans_quantity') +
facet_wrap(~var_type, scales = "free", ncol = 1) + #create one plot for organisms and one for substrate
geom_line()+
labs(color = "Organism/Substrate")
ggsave("oxic_anoxic 2.png")
#fourth switch
name_4='anoxic_oxic 2'
trans4 <- data.frame(out2_t %>% filter(time>340000 & time<348500))
trans4_f <- data.frame(out2_f %>% filter(time>340000 & time<348500))
n_43 <- round(length(trans4$time)*x3/100)
t_43 <- length(trans4$time) - n_43 + 1
n_45 <- round(length(trans4$time)*x5/100)
t_45 <- length(trans4$time) - n_45 + 1
out2_df_2_N_CB <- data.frame(trans4$N_CB)
colnames(out2_df_2_N_CB) <- c("N_CB")
out2_df_2_N_PB <- data.frame(trans4$N_PB)
colnames(out2_df_2_N_PB) <- c("N_PB")
out2_df_2_N_SB <- data.frame(trans4$N_SB)
colnames(out2_df_2_N_SB) <- c("N_SB")
ggplot(aes(x=time, y=trans_quantity, col=species), data=trans4_f) +
ggtitle("Anoxic_oxic 2")+
xlab('Time') +
ylab('Trans_quantity') +
facet_wrap(~var_type, scales = "free", ncol = 1) + #create one plot for organisms and one for substrate
geom_line()+
labs(color = "Organism/Substrate")
ggsave("anoxic_oxic 2.png")
Rolling window analysis
#analysis for 1st
rollingwindow_analysis(trans1, 2, n_15, t_15, name_1)
rollingwindow_analysis(trans1, 3, n_15, t_15, name_1)
rollingwindow_analysis(trans1, 4, n_15, t_15, name_1)
#analysis for 2nd
rollingwindow_analysis(trans2, 2, n_25, t_25, name_2)
rollingwindow_analysis(trans2, 3, n_25, t_25, name_2)
rollingwindow_analysis(trans2, 4, n_25, t_25, name_2)
#analysis for 3rd
rollingwindow_analysis(trans3, 2, n_35, t_35, name_3)
rollingwindow_analysis(trans3, 3, n_35, t_35, name_3)
rollingwindow_analysis(trans3, 4, n_35, t_35, name_3)
#analysis for 4th
rollingwindow_analysis(trans4, 2, n_45, t_45, name_4)
rollingwindow_analysis(trans4, 3, n_45, t_45, name_4)
rollingwindow_analysis(trans4, 4, n_45, t_45, name_4)
Metric based indicators from the earlywarnings package
#analysis for 1st
earlywarnings_plot(out_df_1_N_CB, name_1, x3)
earlywarnings_plot(out_df_1_N_PB, name_1, x3)
earlywarnings_plot(out_df_1_N_SB, name_1, x3)
#analysis for 2nd
earlywarnings_plot(out_df_2_N_CB, name_2, x3)
earlywarnings_plot(out_df_2_N_PB, name_2, x3)
earlywarnings_plot(out_df_2_N_SB, name_2, x3)
#analysis for 3rd
earlywarnings_plot(out2_df_1_N_CB, name_3, x3)
earlywarnings_plot(out2_df_1_N_PB, name_3, x3)
earlywarnings_plot(out2_df_1_N_SB, name_3, x3)
#analysis for 4th
earlywarnings_plot(out2_df_2_N_CB, name_4, x3)
earlywarnings_plot(out2_df_2_N_PB, name_4, x3)
earlywarnings_plot(out2_df_2_N_SB, name_4, x3)
Conditional Heteroskedasticity
#analysis for 1st
ch <- condhet_analysis(out_df_1_N_CB, name_1)
ch <- condhet_analysis(out_df_1_N_PB, name_1)
ch <- condhet_analysis(out_df_1_N_SB, name_1)
#analysis for 2nd
ch2 <- condhet_analysis(out_df_2_N_CB, name_2)
ch2 <- condhet_analysis(out_df_2_N_PB, name_2)
ch2 <- condhet_analysis(out_df_2_N_SB, name_2)
#analysis for 3rd
ch3 <- condhet_analysis(out2_df_1_N_CB, name_3)
ch3 <- condhet_analysis(out2_df_1_N_PB, name_3)
ch3 <- condhet_analysis(out2_df_1_N_SB, name_3)
#analysis for 4th
ch4 <- condhet_analysis(out2_df_2_N_CB,name_4)
ch4 <- condhet_analysis(out2_df_2_N_PB, name_4)
ch4 <- condhet_analysis(out2_df_2_N_SB, name_4)
Detrending fluctuation analysis
#analysis for 1st
dfa_analysis(trans1, 2, n_15, t_15, name_1)
dfa_analysis(trans1, 3, n_15, t_15, name_1)
dfa_analysis(trans1, 4, n_15, t_15, name_1)
#analysis for 2nd
dfa_analysis(trans2, 2, n_25, t_25, name_2)
dfa_analysis(trans2, 3, n_25, t_25, name_2)
dfa_analysis(trans2, 4, n_25, t_25, name_2)
#analysis for 3rd
dfa_analysis(trans3, 2, n_35, t_35, name_3)
dfa_analysis(trans3, 3, n_35, t_35, name_3)
dfa_analysis(trans3, 4, n_35, t_35, name_3)
#analysis for 4th
dfa_analysis(trans4, 2, n_45, t_45, name_4)
dfa_analysis(trans4, 3, n_45, t_45, name_4)
dfa_analysis(trans4, 4, n_45, t_45, name_4)
Drift Diffusion Jump Nonparametrics Early Warning Signals:
#shift 1
ddj11<- ddj_model(out_df_1_N_CB)
plot_ddj_1(out_df_1_N_CB, ddj11, name_1)
plot_ddj_2(out_df_1_N_CB, ddj11, name_1)
ddj12<- ddj_model(out_df_1_N_PB)
plot_ddj_1(out_df_1_N_PB, ddj12, name_1)
plot_ddj_2(out_df_1_N_PB, ddj12, name_1)
ddj13<- ddj_model(out_df_1_N_SB)
plot_ddj_1(out_df_1_N_SB, ddj13, name_1)
plot_ddj_2(out_df_1_N_SB, ddj13, name_1)
#shift 2
ddj21<- ddj_model(out_df_2_N_CB)
plot_ddj_1(out_df_2_N_CB, ddj21, name_2)
plot_ddj_2(out_df_2_N_CB, ddj21, name_2)
ddj22<- ddj_model(out_df_2_N_PB)
plot_ddj_1(out_df_2_N_PB, ddj22, name_2)
plot_ddj_2(out_df_2_N_PB, ddj22, name_2)
ddj23<- ddj_model(out_df_2_N_SB)
plot_ddj_1(out_df_2_N_SB, ddj23, name_2)
plot_ddj_2(out_df_2_N_SB, ddj23, name_2)
#shift 3
ddj31<- ddj_model(out2_df_1_N_CB)
plot_ddj_1(out2_df_1_N_CB, ddj31, name_3)
plot_ddj_2(out2_df_1_N_CB, ddj31, name_3)
ddj32<- ddj_model(out2_df_1_N_PB)
plot_ddj_1(out2_df_1_N_PB, ddj32, name_3)
plot_ddj_2(out2_df_1_N_PB, ddj32, name_3)
ddj33<- ddj_model(out2_df_1_N_SB)
plot_ddj_1(out2_df_1_N_SB, ddj33, name_3)
plot_ddj_2(out2_df_1_N_SB, ddj33, name_3)
#shift 4
ddj41<- ddj_model(out2_df_2_N_CB)
plot_ddj_1(out2_df_2_N_CB, ddj41, name_4)
plot_ddj_2(out2_df_2_N_CB, ddj41, name_4)
ddj42<- ddj_model(out2_df_2_N_PB)
plot_ddj_1(out2_df_2_N_PB, ddj42, name_4)
plot_ddj_2(out2_df_2_N_PB, ddj42, name_4)
ddj43<- ddj_model(out2_df_2_N_SB)
plot_ddj_1(out2_df_2_N_SB, ddj43, name_4)
plot_ddj_2(out2_df_2_N_SB, ddj43, name_4)